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1 Introduction 

rG Feynman Path Integrals ( [8] ) offer a more intuitive alternative to the Schrodinger 

approach in Quantum Mechanics. Green's functions for the Schrodinger opera- 
tor are computed as the limit of sums of exponents of the classical action along 
all possible trajectories connecting two configurations on a lattice, as the lattice 
spacing tends to zero. Path integrals naturally reveal the classical trajectories 
as the limit of "trajectory beams" that contribute the most to the conditional 
probability. This naturally leads to asymptotics that use only some of the 
trajectories - e.g., those "sufficiently close" to the classic trajectories - in the 
J> , integral evaluation. 

However, with the exception of a few special cases of simple Hamiltoni- 
anq^J, Feynman path integrals are notorious for their computational complexity, 
even for reduced "beams" of trajectories. Computing path integrals using the 
classic trajectories only, on the other hand, results in wrong Green's function 
amplitudes (although may still be useful for simple qualitative analysis of quan- 

~f\ turn phenomena as in the double-slit experiment [3]) However, such simplified 

Green's functions can be corrected by applying a normalizing "pre-factor" that 
can be shown to be the regularized operator determinant for the equation in 
variations. Therefore, our ability to compute functional determinants is key 
to successful application of path integrals to semi-classical approximation. A 
natural question to ask in this context is why not use alternative semi-classical 
approximation techniques, such as WKB? ([E]) Methods based on asymptotic 
solution of the Schrodinger equation represent the wave functions (or eigenmodcs 
and energy levels if eigenvalue problem is being solved) as series of powers of h. 
Specifically, the probability density is represented as 

ip(x,t) = A(x)exp ' 



ft 

where 

A = A a + hAi + h 2 A 2 + 



and 



S = Sq + hS\ + h 5*2 



However, it can be shown that for many cases of interest this asymptotic 
representation causes "caustics" - i.e., multi- valued folding singularities - where 
the classical trajectories intersect ([H]). An elegant solution to the problem of 
caustics is delivered by Maslov's Canonical operator ([I2],[H]) that provides 
global-time asymptotic solution by "splicing" together local asymptotics in x 
and ^-representations. The Canonical Operator effectively provides asymptotic 
solution on the phase-space as opposed to the traditional WKB that constructs 
asymptotics in x, thus on the configuration space. The disadvantage of the 
Canonical Operator is that, very much like the traditional WKB, it is essentially 
a device for avoiding caustics, and the splicing procedure - while useful theoret- 
ically - is problematic in practical applications of both Quantum Mechanics and 
wave propagation where short-wavelength approximations are routinely used. 

Functional Integration, on the other hand, is an inherently "global" proce- 
dure that automatically takes care of "widening" or "narrowing" of the con- 
tributing beam of trajectories if all the possible trajectories are included. In 
this work we will review analytical and computational techniques for evaluat- 
ing functional determinants for one and multi- dimensional Hamiltonians. We 
will see that for the one-dimensional and radially-symmetric cases, operator 
determinants can be computed numerically using properties of the Fredholm 
determinants of Storm-Liouville operators f |15]V We provide some numerical 
examples of computing operator determinants and cite references to applications 
in e.g. Quantum Field Theory. However, the case of arbitrary multi-dimensional 
Hamiltonians appears to be significantly more challenging, and apparently no 
computational technique exists that does require the explicit evaluation of at 
least some operator eigenvalues. We offer some tentative ideas for address- 
ing the problem of multi-dimensional Hamiltonians, and discuss the associated 
computational challenges. 

2 Feynman Path Integral 

Solution to the Cauchy problem for the wave function describing motion of a 
single particle of mass to in a potential field V(x) 

mf t = -#.w + vw, (i) 

with the initial condition 

ip(x,t) = ipt (x), t = t 

is given by 

tp(x,t) = / K(x,t;y,t )*l>(y,to)dy (2) 

R" 



where the integration kernel, Green's function or propagator 

i(t-t?) 



K(x,t;y,t') = exp 



-H 



(3) 



has the meaning of the conditional probability of finding the particle at a point 
x at time t given that it was at the point y at time t', where H is the Hamilto- 
nian from the right-hand side of ([lj. The operator exponent in ([3]) is difficult 
to compute for an arbitrary V(x), as it requires the knowledge of eigenfunctions 
and energies. Evaluating the exponent of the Laplacian and the potential sepa- 
rately is straightforward, but exponent of the sum of non-commuting operators 
is not equal to the product of individual exponents. However, commutator of 
two "infinitesimal" operators is an infinitesimal of a higher order, so we can use 
the following approximation 



iAt 



exp- 



P 

2m 



V(x) 



exp 



iAt p 2 
h 2m 



exp 



iAt 



V(x) 



(4) 



for the propagator over a small time period At. Using the semigroup property 
of propagators and substituting Q in © with tp(x,to) — S(x), we obtain, after 
passing At — > 0: 



K(x,t;y,t') 



lim 

JV->oc 
JV-1 

n 



exp 



R"-l 



dpj dxj 
2nh 



. N-l 
3=0 



H{p 3l x )At 



dpp , 
2nh 



(5) 



where H is the symbol of the Hamiltonian and we have substituted the free- 
particle kernel due to the first exponential in the right-hand side of ((4]). Note 
that ([5]) is the limit of summation over trajectories in the phase space (p, x) and 
the integrand is the Legandre transform of the Hamiltonian into the Lagrangian 
(PQ) times the time step. The equivalent configuration space expression is 



K(x,t:y,t') = lim % 2 , 

n^oo V 2nihAt 



exp 



R^- 1 



■ N-l / / N 2 

I x— v I m ( Xj+i — Xj 



E 

3=0 



At 



V( Xj 



At 



dpp 
2nh 



N-l 



n dx 3 



i=i 



(6) 



where the Fresnel integral formula was used for p-integrals. Note that if imag- 
inary time is used, the procedure remains the same but the Fresnel integral is 
replaced with a Gaussian one. Formula (|6]) is interpreted as the summation 
of the exponentiated classical action ([T]) times j- over arbitrary discretized 
trajectories x(tj) = Xj. 



3 Semi-classical Green's Function Approxima- 
tion 

Since we are interested in semi-classical approximations and would like to limit 
the summation in ([B]) to a narrow band of trajectories around the classical tra- 
jectory, let us estimate the variation of the classical action due to perturbations 
of the classical trajectory. Assuming that the perturbation is zero at t and t' , 
replacing the sum in ([6]) with an integral from t' to t and using integration by 
parts we get: 

/ (f \A> - V(,)) « + / (f I A 2 f - igM^A^) „ + (|A,r), 

V t' 

(7) 

where z(t) is the classical trajectory and Az is the perturbation. Substituting 
(J7J) in ([5]) we see that i/ie /irst order correction to Green's function computed 
over the classic trajectories is given by the multiplicative factor 



lim / exp 

JV-i-oo J 

R N-1 



_ / f n/U-.T-'- f^ Az'Az*),// 



2h I \ dxidx k 



N-l 

Y[ dAzj. (8) 



The exponent in ([5J is a quadratic form corresponding to the Hessian of the 
Harmonic oscillator ("Jacobi" matrix). Integral <j8j) is similar to the Fresnel 
integral (or a Gaussian integral if time is changed to imaginary) and is equal to 

,«, (9) 

where |e| = 1 and 

d 2 d 2 V 

J= - m ^--^W- (10) 

Note that operator Hessian (fTO)) maps vector-functions but in the ID case this 
is a simple Storm-Liouville operator 

J=-m*.-V»(z). (11) 

4 Regularized Operator Determinants and Op- 
erator Zeta Function 

Formula (|9]) was not proved but introduced by analogy with finite-dimensional 
quadratic forms. For arbitrary differential operators (|10[) , (JTTJ) we will have to 



define operator determinant that would match the conventional determinant in 
the finite-dimensional case. One such definition utilizes Operator Zeta-function 

my- 

+00 
0(r) = £ A « T < (12) 

where X n n = 1,2,... is the discrete spectrum of J. Given (|12p . operator deter- 
minant is defined as 

dctJ = cxp-Cj(O). (13) 



Note that definition (|13|) coincides with the finite-dimensional determinant if 
J has a finite discrete spectrum and the summation in (|12[) is limited to non- 
zero eigenvalues. The Hessian operator for typical potentials may not have 
a simple discrete spectrum unless some boundary conditions are imposed on a 
finite domain (e.g., Dirichlet boundary conditions corresponding to the infinitely 
high potential barriers at the boundary of interest ([H]). We assume that at 
least for computational purposes such boundary conditions have been imposed, 
effectively cutting off far-field potentials. 

If operator (TTOj) eigenvalues are known, then formula (TT3l) can be used to 
compute the determinant after the operator zeta function is analytically contin- 
ued to zero. The last step is necessary because (J12I) does not define £j at zero. 
The computation can be performed using contour integration ([lip. However, 
the most difficult part - computation of the operator spectrum - is obviously 
not addressed by this procedure. 



5 Computation for one-dimensional and Radially- 
symmetric Hamiltonians 

In the simplest ID case, operator determinant for (TTTJ) can be computed using 
Gelfand-Yaglom theorem ([9], [7], [4]): 



If Dirichlet boundary conditions are 
terminant of the operator (fTTj) actinj 


imposed on an interval 
I on the corresponding 


[a, b], then de- 
Sobolev space 


lem 


given by the y(b), where % 
Jy = 0, 


(x) is 

y(a) 


the solutior 
= 0, y'(a) = 


1 to the initial-value prob- 
= 1. (14) 



Problem ([14]) can be easily solved using e.g. the Runge-Kutta numerical 
method ([2]) for arbitrary V"(x). As a demonstration of this method, we will 
compute the values of the Riemann's zeta function ([16 J at integer points 2 and 



41- Using the formula ([7]) 

fe=i 

we compute the left-hand side for a set of values of A, assuming that Dirichlet 
boundary conditions are imposed on [0,7r] and letting the potential be equal 
to zero. Eigenvalues of J are then trivial to compute and are equal to k 2 , k = 
1,2,.. ., hence the operator zeta function is simply 0( T ) = C(2t). Now fitting 
the tabulated values of the left hand side of (|T5|) with a 4th order polynomial, we 
get 0(1) — C(2) ~ 1-6 and (j(2) = £(4) w 1.1 - both results are in reasonable 
agreement with the exact values (16 ). Of course, the accuracy of this algorithm 
is limited only by the number of points used in the fittingfl, while the accuracy 
of the operator determinant evaluation is limited only by the Runge-Kutta time 
step. 

The Gelfand-Yaglom proposition, although not applicable directly to multi- 
dimensional operators, is part of a similar technique developed for radially sym- 
metric operators ([5]). An application of the radially-symmetric operator deter- 
minants to fluctuation determinant for false vacuum decay is presented in [1]. 
Note, however, that the same problem is solved in [B] using WKB. 

6 Way forward for Multi-dimensional Hamilto- 
nians 



For non-symmetric multi-dimensional Hessians (|10[) it is natural to ask if the 
operator determinant can be computed using some iterative updating procedure 
based on e.g. one-dimensional operator projections. For example, in case of a 
two-dimensional potential pit, solving an initial-value problem similar to (|14p 
along classic trajectories connecting random boundary points would yield the 
operator determinants of one-dimensional operator projections. While these 
one-dimensional determinants are obviously connected with the determinant of 
J, it is unclear how to reconstruct the latter from the former. 
A more traditional approach is based on 

• estimating operator eigenvalue asymptotics; 

• estimating, from the analytical extension of ([i"2j) . the error in (fT3|) due to 
the discarded highest eigenvalues^, and estimating the required number of 
the smallest eigenvalues required to achieve the desired accuracy; 

• computing the required number of eigenvalues using an iterative method, 
e.g. based on Lancsoz iterations ([TO]). 



2 value at 1 is infinity 

3 we used just 10 values with built-in Matlab polynomial fitting 

4 note that due to the analytical extension to zero and differentiation, this error may be 
larger than that of l|12[ l 



Note that for discretization grid sizes of sw 10 3 in each dimension, the corre- 
sponding sparse numerical matrices for e.g. non-symmetric 3D potentials have 
dimensions of 10 9 x 10 9 . Unless only a few initial eigenvalues are required and a 
quick convergence can be expected, this method may be impractical, especially 
in comparison with traditional perturbation and asymptotic techniques. 
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